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ABSTRACT 


DESIGN OF I NFRASOUNO-UETECTION SYSTEM 
VIA ADAPTIVE LMSTDE ALGORITHM 

Camille S. Khalaf 
Old Dominion University 
Director: John W. Stoughton 

A proposed solution to an aviation safety problem is based on pas- 
sive detection of turbulent weather phenomena through their infrasonic 
emission. This thesis describes a system design that is adequate for 
detection and bearing evaluation of infrasounds. An array of four sen- 
sors, with the appropriate hardware, is used for the detection part. 
Bearing evaluation is based on estimates of time delays between sensor 
outputs. The generalized cross correlation (GCC), as the conventional 
time-delay estimation (TDE) method, is first reviewed. An adaptive TDE 
approach, using the least mean square (LMS) algorithm, is then discuss- 
ed. A comparison between the two techniques is made and the advantages 
of the adaptive approach are listed. The behavior of the GCC, as a Roth 
processor, is examined for the anticipated signals. It is shown that 
the Roth processor has the desired effect of sharpening the peak of the 
correlation function. It is also shown that the LMSTDE technique is an 
equivalent implementation of the Roth processor in the time domain. A 
LMSTDE lead-lag model, with a variable stability coefficient and a con- 
vergence criterion, is designed. This model is employed in an automatic 
scheme developed for the sensor array. The software and hardware system 
parameters are derived and determined. The effectiveness of the system 
is illustrated through simulation and field testing. 
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CHAPTER 1 


INTRODUCTION 
1.1 Motivation 

Current research programs have lead to significant advances in 
ground-based and airborne equipment for providing information relative 
to severe turbulent weather. However, most of these programs are still 
in their experimental stages and very few are operational. There con- 
tinues to be a serious aviation safety problem associated with aircraft 
operations in the vicinity of severe storms. 

In 1981, general aviation aircraft numbered more than 200,000 and 
flew more than 40 million hours [1]. General aviation operations 
accounted for 662 fatal accidents from all causes, with 1,265 fatalities 
(FAA, 1981). Informal accident cause/factor statistics from the Nation- 
al Transportation Safety Board for 1981 indicate that weather caused, or 
was a factor, in 40 percent (289 cases) of the U. S. general aviation 
accidents. Earlier statistics indicated that turbulence is the largest 
single cause of weather-related air carrier accidents in the U. S. From 
1962 to 1974, turbulence was either a cause of or a contributing factor 
in 189 of 450 weather-rel ated cases [2], 

Tne two types of turbulence usually encountered are clear-air tur- 
bulence (CAT) and thunderstorm-related turbulence. CAT, a problem for 
all aircraft, cannot be seen because it usually has no cloud signature. 
It may develop in a standing wave caused by air moving over mountainous 
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terrain, and is frequently associated with shear-induced Kelvin- 
Helmholtz atmospheric waves occurring in a statistically stable atmos- 
phere [3]. Accidents caused by CAT are not as serious as the ones 
related to turbulence associated with thunderstorms. A CAT accident may 
result in discomfort, injuries, aircraft damage, and/or unscheduled 
landing. Thunderstorms and other convective clouds are critically 
important sources of low-altitude turbulence and wind variability. Many 
produce strong downdrafts that transport air downward, which then 
spreads out rapidly over the ground. This mechanism, if encountered 
during take off or landing of aircraft, may result in serious if not 
fatal accidents. 

One of the pressing aviation safety problems is that of providing 
the pilot with information needed to avoid turbulence hazards which 
exceed the design capabilities of the airplane. A proposed solution to 
this proDlem through passive detection of turbulence is presented in the 
next section. 

1.2 Turbulence Detection 

The primary technique for detecting turbulence and storm cells is 
Doppler radar. This technique requires trie presence of reflective 
particles such as precipitation or dust for meterological applications. 
Doppler radar proved ineffective in cases of CAT or developing storm 
cells due to the absence of reflective particles in these phenomena. 

The proposed technique consists of passive detection of large-scale 
patches of turbulence in the Earth's atmosphere through infrasonic 
emissions. The infrasonic technique offers the advantages of being (1) 
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passive, (2) inexpensive, and (3) inherently more sensitive to atmos- 
pheric disturbances than electromagnetic or optical techniques. The 
assumed utilization of the infrasonic technology is illustrated in Fig- 
ure 1.1. Infrasonic emissions from the patch of turbulence are detected 
at two stations, each containing an array of four sensors. By means of 
real-time signal processing, the direction of the turbulent source, 
based on time delays between sensors, is determined at each station. 

The turbulence is then located by tri angul ation between the two sta- 
tions, This information is relayed from a control center to a flight 
services advisory, and from there to the pilots of approaching aircraft. 
A series of such stations would provide early warning along the length 
of domestic air traffic routes. 

It is informative at this point to mention that the first step 
toward realization of the infrasonic technique was the development of a 
unified acquisition system for acoustic data [4] by Dr. Allan J. 
Zuckerwar and Mr. Harlan K. Holmes at NASA Langley Research Center. 

This system served an important role in understanding the detection 
problem during the early stage of this research. 

The del ay-to- angle conversion needed in this passive approach is 
illustrated in Figure 1.2. If is the time delay between the 
arrival of the infrasonic signal at two spatially separated sensors, 
then the direction of the turbulent source is given by 

e = cos~ I ( ) (1.1) 

z 
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TURBULENCE 



Figure 1.1. Early warning system. 
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where V is the sound velocity (1080 ft/sec at 0°C), and z is the 
distance between the two sensors. For a given delay, the solution of 
equation (1.1) can be interpreted in two ways as shown in the figure. 

The signal associated with the false direction is usually referred to as 
the "ghost" signal. The presence of the ghost-signal requires a minimum 
of three sensors in the array to uniquely determine the source direc- 
tion. 

1.3 Research Objectives 

The infrasonic technique, proposed to provide pilots with an early 
warning of turbulence, was described in the previous section. The three 
aspects that are essential to its success are the following: 

1. emission and propagation of infrasounds from turbulent weather 
phenomena; 

2. detection of infrasounds and bearing evaluation of their 
sources at every station; 

3. locating the sources by triangul ation between stations via a 
communication network that is also responsible for transferring 
the information to pilots. 

The first aspect is out of the designer's control and is based on the 
theoretical predictions of infrasonic emission and propagation. Tne 
third aspect constitutes the commercial phase of the program once the 
performance of the individual stations proves satisfactory. 

The goal of this research is the overall system design that is 
adequate for implementing tne second aspect of the infrasonic technique. 
In this context, the system should function as an infrasonic detector. 
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time-delay estimator, and delay-to-angle converter. The detection part 
involves the sensor array and corresponding hardware blocks while the 
delay-to-angle conversion is accomplished by a simple software routine. 
The main design issue is that of an appropriate signal processing tech- 
nique for time-delay estimation. The GCC method and the LMS parametric 
technique will both be reviewed. More specifically, the GCC method 
using a Roth weighting function, and the adaptive implementation of the 
LMS technique, will be examined. Design decisions, regarding the TDE 
technique as well as the hardware system, will then be made according to 
the anticipated signal characteristics. 

In Chapter 2, the theory behind TDE techniques is reviewed. Based 
on this theory and the anticipated infrasounds, the system design is 
presented in Chapter 3. Chapter 4 evaluates the system through simu- 
lation of the basic TDE algorithm response and discussion of the system 
performance in the actual field. Conclusions and future research are 
presented in Chapter 5. 



CHAPTER 2 


THEORY 

2.1 Introduction 

The problem identified in chapter one is a fundamental passive 
sonar signal processing problem in which delays between the times of 
arrival of the pertinent acoustic waves at four sensors are to oe esti- 
mated. This chapter will discuss the theory behind TDE techniques. 
Section two reviews the conventional generalized cross correlation 
approach (GCC) while section three reviews the Roth processor specifi- 
cally. Of particular interest is a parametric approach through widrow's 
adaptive filter. The filter structure is reviewed in section four while 
its application to TOE is presented and thoroughly investigated in the 
fifth section. The chapter is concluded in section six by presenting 
the advantages of the adaptive least mean squared time delay estimation 
(LMSTDE) method. 

2.2 Generalized Cross Correlation Approach 
The GCC approach to time delay estimation has been discussed by 
many investigators. Well known references include papers written by 
Knapp and Carter [5], and Hassab and Boucher [6]. This section only 
reviews this approach through Figure 2.1. x x (t) and x 2 (t) are sam- 
pled at two spatially separated sensors and then fed to a basic cross 
correlator. The Dasic cross correlator, as discussed by Papoulis [7], 
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Figure 2.1. Generalized cross correlation. 
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computes the cross correlation function, R (x), between Xi (t) and 

xi x 2 

xz(t) by means of the inverse Fourier transform of their cross pow er 
spectrum, G (f ) . Then the delay estimate is simply the abscissa 

Xi X2 

value at which the cross correlation function peaks. In order to im- 
prove the accuracy of the delay estimate, a linear filter, w(t), is 
convolved with the output of the cross correlator. A peak detector 
routine is then used to determine the abscissa value of the peak in the 
filtered correlation function. In practice, a frequency weighting W(f) 

= F [w(t)], where F[«] denotes the Fourier transform of [*J, which is 
equivalent to w(t), is applied to the cross power spectrum prior to 
taking the inverse Fourier transform. This frequency weighting replaces 
the linear filter so that all computations, except for the peak detector 
routine, are done in the frequency domain. A discussion of this 
weighting, after mathematically model ing the system, is presented next. 

A signal emanating from a remote source and monitored in the pres- 
ence of noise at two spatially separated sensors can be mathematically 
modeled as 


xi (t) = Si (t) + mi(t) (2.1a) 

X 2 (t) = aSi(t + x d ) + m2(t) (2 . Id) 

where Si(t), mi (t) and m 2 (t) are real, zero-mean, jointly stationary 
random processes, x d denotes time delay and a an attenuation factor. 
Signal S x (t) is assumed to be uncorrelated with 1 % (t) and m 2 (t). 
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The cross correlation between x 1 _(t) and x 2 (t) is related to the 

cross power spectral density function, G (f), by the Fourier trans- 

X 1 x 2 

form relationship 


R 

XI X2 


(t) = 


00 

/ G (f) e J ^ fT df 

—oo XI X2 


( 2 . 2 ) 


In practice, only an estimate G* (f) of (2.2) can be obtained from 

x i x 2 

finite observations of xi(t) and X 2 (t). Consequently, the output of 
the basic cross correlator, when no weighting is used, is: 


oo 

R* (t) = / G’ (f) e j2irfT df (2.3) 

x l x 2 ioo X 1 x 2 


It is informative at this point to examine the shape of the cross 

correlation function R (t), so that the necessity of using a 

X 1 x 2 

weighting function can be justified. An expression for R XlX2 ^ T ^ can 
be obtained from (2.1), using the expectation operator E[* ] , as 

V 2 (t) = E[xi(t) Xz(t ' T] = ° R SiS 1 ( T - T d) + S"2 (T) (2 - 4) 

The Fourier transform of (2.4) gives the cross power spectrum 


G 

XI X2 


(f) 


aG c c (f) 
Si Si 


— j 2 it f t 


+ G (f) 
mi m 2 


(2.5) 
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If m 1 (t) and m^t) are uncorrelated, G mi ni 2 ^) = and s '’ nce 
multiplication in one domain is a convolution in the transformed domain, 
it follows from (2.5) that 

R (t) = aR - (x) 8 6 (t-t ) (2.6) 

x 2 Si Si d 

-j2irfx d 

wnere ® denotes convolution and 5(t-t^) = F- 1 [e ]. one 

interpretation of (2.6) is that the delta function has been spread or 
"smeared" by the Fourier transform of the signal spectrum. If Sj(t) 
is a white noise source, then its autocorrelation function is a delta 
function and no spreading takes place. However, for most practical 
applications this is not the case and spreading acts to broaden the peak 
of the cross correlation. 

To minimize the spreading effect, many weighting functions have 
been proposed in the literature (see Table 2.1) to operate on the cross 
power spectrum given in (2.5). With a general weighting, Wg(f), the 
estimate of the generalized cross correlation becomes 

oo 

R 9 (t) = / Wg(f) G (f) e j2ufT df (2.7) 

x l x 2 _oo x l x 2 

Wg(f) should be cnosen to ensure a large sharp peak in R 9 (t) rather 

x l x 2 

than a broad one in order to obtain a good time-delay resolution. 

2.3 Roth Processor 

The weighting functions found in the literature are listed in Table 
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2.1, where the notation y(f) has been used for the coherence function 
which is defined as 


Y12 


(f) = 


G 2 (f) 

v v ' ' 


x l x 2 


(f) G v „ (f) 


Xi X 


1*1 


x 2 x 2 


The selection of Wg(f) ..to optimize certain performance criteria has 
been studied by several investigators (see, for example [8]). The pur- 
pose of considering the Roth processor is that it is equivalent to the 
time domain weighting involved in the adaptive LMSTDE approach. This 
equivalence will be shown in section five and will be helpful in gaining 
insight to the adaptive method. 

The weighting proposed by Roth [9] is 


W f > - 


Substituting for Wg(f) in (2.7) yields 


R v v ( T ) * J 


XiX 2 


(Wiener-Hopf) filter [10] 


H(f ) = 


v 

(f) 

1 

yields 


G x x 

x i x 2 

(f) 

S*, 

(f) 

se res 

ponse 

G x X 

S X 1 X 2 

(f) 

G 

x i x l 

(f) 


>32^ f x 


df 


( 2 . 8 ) 
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Processor Name 


Weight 



1 Smoothed coherence transform 
2 Phase transform 
3 Maximum likelihood 


Reference 

[ 8 ], [ 7 ] 

[ 9 ] 

[ 10 ] 
[ 10 ] 
[ 5 ] 

[ 5 ] 
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which "best" approximates the mapping of Xj (t) to x 2 (t). If m 1 (t) * 
0, then G v v (f) = G<- <• (f) + G (f), and equation (2.6) becomes for 

Xj A1 ^1^1 l"l 

the Roth processor 


R (t ) = 5 (t-t ) b / 
XI X2 d 


aG (f) 
5 1 5 1 


G, c ( f ) + G m m 

5i 5i mi mi 


e j2irfT df(2. 9) 


From (2.9) we can conclude that only when G^ ^ (f) is negligible or 
when it equals any constant times G^ (f), the spreading does not 


occur and R (t) becomes a delta function. However, as can be seen 
x l x 2 

from the integral part in (2.9), the Roth processor has the desirable 

effect of suppressing those frequency regions where G m m (f) is large 

and G (f) is more likely to be in error. 
x l x 2 

The generalized cross correlation approach suffers regardless of 
which weighting function is used, from two basic facts. First, it 
relies on a sequence of fast Fourier transform (FFT) computations that 
tends to be time consuming. Second, as seen in Table 2.1, it requires a 
priori knowledge of the signal and noise statistics to implement the 
specific weighting. In passive detection problems, this information is 
unknown and if it were to be estimated, it would increase the complexity 
of the process and the time involved. The next section presents the 
adaptive LMSTDE algorithm, which is equivalent to the Roth processor, 
and is able to overcome both difficulties found in the generalized cross 
correlation method. 
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2.4 Adaptive Filters 

This section will review the basic structure of adaptive filters 
discussed by Widrow [11], while section 2.5 presents the application of 
the filter in TDE problems. 

A signal filtering approach using an adaptive filter is in some 
sense a self-designing (really self-optimizing) process. The adaptive 
filter described here bases its own "design" (its internal adjustment 
settings) upon estimated or measured statistical characteristics of the 
input and output signals. The statistics are not measured explicitly; 
rather, the filter design is accomplished in a single process by a re- 
cursive LMS algorithm that automatically updates the system coefficients 
with the arrival of each set of data samples. Figure 2.2 illustrates 
schematically the adaptive filter used in this case as a linear combina- 
torial system. The filter consists of a set of variable weights (filter 
coefficients) whose input are the sampled input signals, a summer to add 
the weighted signals, and an algorithm to adjust the weights automatic- 
ally. The impulse response of such a discrete system is completely 
controlled by the weight settings. The adaptation process automatically 
seeks an optimal filter impulse response by adjusting the weights using 
gradient techniques to minimize the mean-square-error function. 

The analysis of the adaptive filter can be developed by assuming 
that the input signals are statistically stationary random processes. 

Let the nth set of input signals be a vector X(n) of length N, 

j 

X (n) = [x x (n) x 2 (n) ... x N (n)] 
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Figure 2.2. Adaptive linear combinatorial system. 
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where X denotes the transpose of X. Let the set of weights, at the 

nth time, be designated by the vector w"'"(n) = [w x (n) w 2 (n) ... w N (n)]. 

N 

The nth output signal is y(n) = l w,(n) x^n). This can be written 

i=l 1 1 

in matrix form as 


y(n) = W T (n) X(n) = X T (n) W(n) . 


( 2 . 10 ) 


Denoting the desired response by d(n), the error at tne nth time is 


e(n) = d(n) - y(n) = d(n) - W (n) X(n). 


( 2 . 11 ) 


The square of tne error is 


e 2 (n) = d2(n) - 2d(n) X T (n) W(n) + W T (n) X(n) X T (n) W(n) . (2.12) 


The expected value of e 2 (n) is 


E [s 2 (n) ] = cP(n) - 2 vj x W(n) + W T (n) V xx W(n). (2.13) 


where (• ) denotes the mean value of (*), or E[*]. The vector V dx is 
the cross covariance between d(n) and x(n) and is defined as 


d(n) xi(n) 

V dx = E[d(n) X(n)] = E d(n) X2(n) 

*d(n) x N (n) 


(2.14) 
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and V denotes the auto-covariance matrix of X(n) 

AA 



* E [X(n) X T (n) ] = E 


"x l ( n ) x x (n) 
x 2 (n) xx(n) 


Xi(n) x 2 (n) . . . 
x 2 (n) x 2 (n) .. 


(2.15) 


V n ) x N ( n ) 


It may be observed from (2.13) that for stationary input signals, the 
mean-square error is a second-order function of the weights. Thus, the 
mean-square-error function may be viewed, as suggested by Widrow [11], 
as a "performance surface" for the adaptive process that has a unique 
stationary point (minimum) which can be sought using gradient tech- 
niques. 

The gradient at any point on the performance surface can be obtain- 
ed by differentiating the mean-square-error function (2.13) with respect 
to the weights. The gradient is 


V[e 2 (n)] = -2 V dx + 2 V xx W(n) 


(2.16) 


The "optimal" weight vector, W^, that yields the least-mean-square 
error, is found where the gradient is zero. Accordingly, 


V, = V W, MC 
dx xx IMS 


“LMS^xx)- 1 V <ix 


or 


(2.17) 
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where [• j- 1 denotes the inverse of [•]. Equation (2.17) is the Wiener- 
Hopf equation in matrix form [12]. Solving (2.17) for the optimum 
weight vector (the minimum point on the performance surface) can present 
serious computational problems. However, it will be shown that the 
adaptation process tries to find an exact or an approximate solution to 
the Wiener-Hopf equation by using the LMS algorithm with less computa- 
tional complexity. 

When using the LMS algorithm, changes in the weight vector are made 
along the direction of the estimated gradient vector. Accordingly, 

W(n+1) = W(n) + K s V [ 7 2 (n)] (2.18) 


where 


W(n) A weight vector before adaptation 
W(n+1) 4 weight vector after adaptation 

A scalar constant controlling the rate of convergence and 
stability (K $ < 0) 

V'[e 2 (n)]4 estimate of gradient of E[e 2 ] = e 2 with respect to W, 
with W = W(n) 

One method for obtaining the estimated gradient of the mean-square-error 
function is to take the gradient of a single time sample of the squared 
error. That is, 
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V'ip(n)] = V [s 2 ( n ) ] = 2 e (n) V[e(n)] (2.19) 

From equation (2.11), 

V[e (n) ] = V [d ( n) - W T (n) X(n)] = - X(n) 

Thus, 

V'[7 2 (n)] = - 2 e(n) X(n) (2.20) 

It can be shown that the gradient estimate of (2.20) is unbiased, so 
that 

E { V [7 2 (n)j = V [7 2 (n) ] 

Substituting (2.20) in (2.18) yields 

W(n+1) = W(n) - 2 K g e(n) X(n) (2.21) 

and the next weight vector is obtained by adding to the present weight 
vector the input vector scaled by the value of the error. This is the 
LMS algorithm. 

Next, the expected value E[W(n)] of the weight vector after a 
large number of iterations will be shown to converge to the Wiener solu- 
tion given by (2.17). For this purpose, assume that the time between 
successive iterations of the algorithm is sufficiently long so that the 
input vectors X(n) and X(n+1) are uncorrel ated. Taking the expected 
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value of both sides of (2.21) we obtain a difference equation in the 
expected value of the weight vector: 


E[W(n+l) ] = E[W(n) ] - 2 K $ E { X(n) [d(n) - X T (n) W(n)]} 

= C I+ 2 K s V xx J E[W(n)] - 2 V dx (2.22) 

.* * 

where I is an identity matrix. With an initial weight vector W(0), 
n+1 iterations of (2.22) yield 

, i n . 

E[W(n+l) ] = [1+2 K s V xx ] n W(0) - 2 K $ E [1+2 K $ V^] 1 V dx (2.23) 

Equation (2.23) can be put in diagonal form by using the normal-form 
expansion of the matrix V . That is, 

aA 


V = Q- 1 A Q 
xx H 

where the diagonal matrix of eigenvalues is A, and the square matrix 
of eigenvectors is the modal matrix Q. Equation (2.23) may now be 
expressed as 

n+i n i 

E[W(n+l) = [1+2 K $ Q- 1 AQ] W(0) - 2 K $ [1+2 K $ Q- 1 AQ] V dx 

= q- 1 [1+2 K A] n+1 qw(0) - 2 K q- 1 S [1+2 KA]^ V Hy .(2.24) 
s s i=Q s ox 
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Consider the diagonal matrix [1+2 K S A]. As long as its diagonal terms 
are all of a magnitude less than unity, 

Aim [1+ 2 K s A] n+1 * 0 
n + °° 


and the first term of (2.24) vanishes as the number of iterations i n- 
creases. The summation factor in the second term of (2.24) becomes 

n . 

Aim z [1+2 K A] 1 = - — — A- 1 
n + “ i=Q 2 K s 

where the formula for the sum of a geometric series has been used. That 
is. 


E ( 1 + 2 K X ) 1 = i = -ll — 

i=0 b l-(l+2 K S X) 2 K S A 

Thus, in the limit, equation (2.24) becomes 


Aim E[W(n+l)] 
n + 00 


Q-^- 1 Q 




which is the Wiener-Hopf solution in (2.17). 

Convergence of E[W(n+l)] to (2.17) is obtained if and only if the 
diagonal terms of [1+2 K $ A] all have magnitudes less than unity, and 
since all eigenvalues in A are positive (the auto-covariance matrix, 
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V , is positive definite), the bounds on K are given by 

A A J 


I 1 + 2 K $ Xj < i or — < K $ < 0 

X m 


where X is the maximum eigenvalue of V . The convergence condition 
m xx 

can be related to the total input power as follows: 


X < trace [V ] = Z E[x1p = Total input power 
rn xx .• » 


Therefore; 


-1 


N 

s EDa 

i=l ! 


< K < 0. 

s 


For a slow, precise adaptation is usually chosen such that 


» IK. 


N 

z ECxp 

i=l 1 


It is believed that tne assumption of independent successive input 
vectors used for the convergence proof is overly restrictive. Griffiths 
[13J has shown that adaptation using highly correlated successive sam- 
ples converges to the Wiener solution, but leads to higher steady-state 
mean-square error. Thus, we can conclude that for short-term stationary 
signals, with the feedback estimate, K s , being bounded by (2.25), the 
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weight vector, after a large number of iterations, is expected to con- 
verge to the weiner solution that best maps the input x(n) to the 
desired response d(n). 

2.5 Adaptive LMSTDE Approach 

The LMS adaptive filter described in the previous section has 
been widely applied in situations where the statistics of the inputs are 
either unknown or partially unknown. Some applications include noise 
cancelling [14], line enhancing [14], prediction [11], spectrum analysis 
[13], and adaptive array processing [15]. A recent application of the 
filter is time delay estimation by F. Reed, P. Feintuch and N. Bershad 
[16]. This application is demonstrated in Figure 2.3, where the adap- 
tive filter has a slightly different structure. The filter to be con- 
sidered here has only two inputs: a primary input x(n) and a second- 

ary input (desired response) d(n). The primary input is fed to a tap- 
ped delay line to generate the adaptive filter input signals. In this 
case, the input vector at the nth iteration becomes 

X T (n) = Dd (n) X 2 ^ n ) •• x n (")l 

= [x(n) x(n-l) ... x(n-N+l)] 

(for a filter of length N) 

and the weights are updated with the arrival of each new data sample 
x(n) . 
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Let xi (t) and x2 (t), in Figure 2.3, be given as in equation 
(2.1) with the same assumptions. Let us consider the discrete time 
version of (1) for the following analysis. That is, 

xi (n) = Si (n) + mi(n) 

and X 2 (n) = S 2 '(n) + m 2 (n) = a Si(n-D) + m 2 (n) 

where D is a positive integer less than N, representing the discrete 
time delay between S x (n) and S 2 (n). 

The adaptive filter inputs are then given by 

x(n) = xi (n) = Si(n) + mi(n) (2.26a) 

and d(n) = X 2 (n) = a Si(n-D) + m 2 (n) (2.26b) 

The adaptive filter, with inputs given by (2.26), can be thought of as a 
system attempting to insert a delay equal to the propagation delay be- 
tween the two sensors in the primary filter input, x(n), aligning the 
signal component in time prior to subtraction to produce the error sig- 
nal. Hence, one weight of the filter corresponding to the correct 
delay would be unity and all other weights zero. In practice, due to 
tne fact that the filter must interpolate between its discrete taps to 
provide delays that are noninteger multiples of the sample time, the 
weights converge to a shape that is peaked at the correct delay. 



I 



Figure E.3. Modeling TDE using adaptive filter. 
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Therefore, determination of the delay, as in the GCC method, requires 
estimation of the peak of the adaptive filter weight vector. 

To examine the shape of the weight vector, let us consider the 
frequency domain window corresponding to the least-mean-square set of 
weights, W LMS = ( v xx ) _1 V dx* That is > 


G dx (f) 

W LMS^ " F ^ W LMS^ k ^ 


G xx( f > 


For d(n) and x(n) given by (2.26), this becomes 


w LMS (f) 


a G<- (f) 


-J2 7T fD 


+ V, (f) 


(2.27) 


Converting back to time domain and using the convolution theorem, the 
kth weight at the nth iteration is 


W LHS (k) = F- 1 [e' j2 * fD ] » F-![. . ta5 J S . L -— ] 


G s lSl < f ) + S% (f) 


- F -‘ [e - j2 " f0 ] * / [- 


“ Sc c (f) 
■>1 ^1 


G S 1 S 1 + 


J e J ’ 2lTfK df (2.28) 


where 0 < K < N-l. Equation (2.28) is to be compared with equation 

(2.9), where W (K) is a discrete time version of R (t) obtained 
LM5 *1 X 2 

using the Roth processor. It is important at this point to note that 
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equations (2.27) and (2.28) do not depend on G^^ (f ) . Therefore, as 
suggested by Ahmed and Carter [17], it is desirable in this adaptive 
approach to have the signal with lower S/N ratio as the secondary input 
to the filter. From equation (2.27) we can see that when 


m 1 m x 


(f) = 0 


(2.29a) 


or. 


G (f) = const, x G c c (f) (2.29b) 

nil 


we have rt LMS (f) = 0 e~ j2irfD , where 3 is a real constant. In this 
case, equation (2.28) becomes 


W lms (K) = 6 x F-l [e' j2lTfD ] 


(2.30) 


Let us examine equation (2.30) for two types of input signals. 

First, consider the case when the input signals are broad-band, i.e. 
f f 

— -< f < — , where f. denotes the sampling rate. Evaluating (2.30) 
2 2 s 

for this type of signal yields 


LMS 


(K) .*/* e -J* fD * e^df = 0 I"' 2 e-^ f0 x e 
- - f s/2 


j2irfK 


df 


= 3 sine {it [K-D]} 


(2.31) 
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where sine (• ) = s j n i.*J. and 0 < K < N-l. In this case the weight 

. (•) 

vector has the shape of a sine function peaking with amplitude equal to 
one when K-D =0 or K = D. The first zero crossing of (2.31) happens 
for 


n (k-D) = ±tt=> K = 0 ± 1 


Hence, the main lobe of the sine function is only two resolutions wide, 
simulating a delta function. However, as in the Roth processor, when 
the condition (2.29) on G m _ (f) is not met, the main lobe of the sine 

}7ij FFK 

function will be spread out by convolving with the integral part of 


(2.28). 

A second relevant case is one in which the input signals are band- 
limited in frequency. Evaluating equation (2.30) for band-limited sig- 
nals leads to a slightly different form of the weight function as dis- 
cussed by Ahmed and Carter [17]. Consider an ideal band-limited signal 
as follows 


be c (uj ) =£ 0 
°1 °1 

Go <-(<») = 0 


for o>o< jw |<u>i 


elsewhere 


where the conversion to the radian frequency, co, was made for the ease 
of analysis only. For this type of signal, equation (2.30) becomes 
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W lms (K) - B /* e-> D e>* <b > 2 B /' e + j“ (K -°>4» 

-IT O ) 0 

= 46 o> b sine [o» (K-D)] cos [to (K-D)] (2.32) 

where w h = BLLzigL and w r = 

D 2 c 2 

The peak in the weight function still occurs at the correct delay, that 
is, for K-D = 0 or when K=D, but the main lobe is spread out. That 
can be seen by examining the first zero crossing of equation (2.32). 

The first zero in the sine function occurs for 

oj . (K-D) = ±tr-> K = — + D 

w b 

while the first zero in the cosine function occurs for 

<o (K-d) = ± - * K = — + D 

c 2 2 u 

c 


Since & > <u. , the first zero in the cosine function occurs sooner 

c b’ 

than in the sine and yields a main lobe width of 2 x — — = — . So for 

2u) a) 
c c 

band-limited input signals, the spreading in the peak of the weight 

function is inversely proportional to the band width. For a white input 

signal, w = L and the main lobe width is two resolutions as pointed 
c o 


out earlier. 
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Other issues that should be considered in this approach include 
interpolation, dynamic behavior of the filter, and the filter length. 
The time delay obtained from the estimated discrete weight function, 
Wlm$(K), is an integer multiple of the sampling period while in prac- 
tice the actual delay between the input signals may not be. To avoid 
the error introduced by the discrete delay estimate the weight function 
should be interpolated to obtain its continuous counterpart Wy^njK), 
that is, 


P 

W LMs( n * K ) = 1 W MS^ m ’ K ) sinc (2.33) 

m=-p ><J 

Now, the desired noninteger sample of the estimated delay is given by 
the value of n at which W LMS (n,K) is maximum. The study of the 
filter dynamic behavior by P. Feintuch, N. Bershad, and F. Reed [18], 
has shown a potential to track linearly changing time delays by just 
observing the peak weight move through the adaptive filter tapped delay 
line. It was also shown that the location of the peak of the weights 
lags the true time delay by an amount that depends on the delay rate, 
the signal correlation function, and the adaptive filter time response. 
Finally, the filter length should be chosen carefully to insure fast and 
yet accurate convergence. The lower bound on tne number of weights is 
depicted by the maximum propagation delay between sensors. From equa- 
tion (2.32) we can see that as K increases the weight function 
W LMs( K ) a PP r oa ch e s zero. Thus, the filter length can be safely 
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truncated in practice to the lower bound, in order to reduce the compu- 
tational time. 

This section is concluded by a representative example illustrating 
the capability of the adaptive filter to detect the time delay after 
convergence of the weight vector. Let x x (n) be a stationary random 
process with a correlation function given by 

/ 

R(n) = a^ with -1 < a < 1 

Assume that mi(n) = m 2 (n) =0 and X 2 (n) = xi(n-D) where D is a 
positive integer given by 0 < D < N-l. The adaptive filter inputs are 
then. 


X (n) = Xx(n) 

d(n) = xi (n-D) 

The input vector on the tapped-delay line is 

X T (n) = [x ( n ) x (n-l) ••• x(n-N+l)] 

and the cross-covariance vector, , is 

vj x = E [d(n)X T (n)] = E [d(n) x (n) d(n)x(n-l) ••• d(n) x(n-N+l)] 
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X (n) = [x(n) x(n-l) ••• x(n-N+l)] 


and the cross-covariance vector, V^, would be 


,T 


V dx = E Cd(n)X(n) 3 = E [d(n)x(n) d(n)x(n-l) ••• d(n) x (n-N+1) ] 


= [R ( — D) R(1->D) ••• R(0) ••• R(N-2-D) R(N-l-D)] 


r D D-l 

= [a a 

• • • 

• • 2, • • • 

] a* 

• • 

N-2-D 

ci 

The auto-covariance matrix, V xx , of 

x(n) 

is 




R(0) R(l) 

• • 


R(N-l) 

T 


R(l) R(0) 

• • 



V xx - E[X(n)x'(n)] 

1 “ 

• 

• 

• 



• 

• 

• 



R(N-l) 

• • 


R(0) 



‘1 a a 2 

• • 


a N ' r 



ala 

• • 


• 


= 

a 2 



• 



• 

• 

• 



1 a 



N-l 



a 1 



_ a • • 




' 1 

-a 0 0 

• • • 




-a 

1+a 2 -a 

• • 



a"d = 

0 

-a 1+a 2 




I 0 

0 -a 





• 

• 

0 0 

• • 



• 

• 


• 

• • 

• • 

• 0 

-a 

1+a 2 


0 

0 0 

• • 

0 

-a 


a^-i-U] 


1-a 2 


-a 
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W 


LMS 


0 

0 

0 


0 

0 


0 

1 

01 


i. L. 

D un row, i.e. the weight vector peaked at the 
correct delay. 


2.6 Summary 

The generalized cross correlation approach was reviewed with more 
attention paid to the Roth processor. The adaptive filter and its 
application to TDE problems was analyzed in detail. It was shown that 
the adaptive LMSTDE method is equivalent to the Roth processor and that 
the adaptation process using the LMS algorithm, under stated conditions, 
ensures the convergence of the weight vector to the wiener-Hopf solution 
that presents the best mapping between the two inputs of the filter. 

This mapping resulted in a weight vector peaking at the correct time 
delay with a peak resolution that is a function of the frequency band of 
the signals and their statistical characteri sties. The advantages of 
the adaptive LMSTDE over the GCC method can be summarized as follows: 

- simpler implementation, since no FFT computations are involved 

- no prior knowledge of the statistics of the signal and noise are 
necessary 

- ability to track moving sources 
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- less computational time. 

The advantages over the nonadaptive approach (i.e. direct implementation 
of the Wiener-Hopf equation) are 

- simpler implementation, since no matrix operations are involved 

- the ability to track moving sources 

- the filter length can be easily increased in order to decrease 
the bias in the estimate. 



Chapter 3 


SYSTEM DESIGN 
3.1 Introduction 

The theory behind time-delay estimation techniques has been review- 
ed in Chapter 2. This theory can now be used in the design process of 
the system. Based on the anticipated input signals and assumptions that 
can be made about the statistical characteristics, the choice of a suit- 
able TDE technique is made in section 2. Section 3 then presents the 
implementation of this technique and relevant design issues in the 
specific application on hand. The system specifications and the way 
they are related, derived, or determined are the subject of section 4. 
Finally, the design process is summarized in section 5. 

3.2 Design Assumptions 

The choice of the TDE approach and the assumptions on which it is 
based are the subject of this section. Not much study has been dedi- 
cated to the high frequency band (1 to 20 Hz) of infrasounds in the 
past, especially for sounds related to meteorological events. Thus, 
little is known about the statistics of these infrasounds and assump- 
tions have to be made with regard to recent investigator's predictions. 
These assumptions will be listed below with a relevant discussion about 
each of them. 

a. Single Source. More than one meteorological event can take 
place at one time producing a multiple delay between sensor 
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outputs. The potential to separate delays and identify the 
different sources is feasible but it is outside the scope of 
this study. So the design process is restricted to a single 
source of sound waves. 

b. Stationary Signal. The stationary assumption is safely made 
because the data file spans a time window that is only a frac- 
tion of the event's lifetime. The lifetime of a meteorological 
event can be anywhere from a few minutes, in the case of a 
microburst for example, to several hours in the case of a fron- 
tal system [1], [19]. Tnis time frame is to be compared with 
the file time that, as will be seen later, is 20 seconds. 

c. Uncorrel ated Signal and Noise. Infrasounds emitted by a rela- 
tively remote source cannot be expected to correlate with local 
noise present at the microphone array. 

d. Uncorrel ated Noise at the Different Sensors. This assumption 
can be justified by having a long array foot so that local 
noise at one sensor is not swept across the array.- 

e. Frequency Band. The low frequency band (10 -6 to about 1 Hz) of 
the infrasound spectrum is more related to static pressure 
fluctuations that are part of a turbulent motion. Since in 
this project, the intention is to detect acoustic (propogating 
wave) energy and as little as possible of the turbulent pres- 
sure field, the frequency band should be in the higher part of 
the spectrum. The l-to-16-Hz band has been experimentally 
studied by Eric S. Posmetier [20]. He concluded that "the 
l-to-16-Hz infrasound may be radiated by clear air turbulence 
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and may be a basis for a remote passive detection system." 

Based on the physics of infrasounds and Posmetier's prediction, 
a signal frequency band of l-to-20-Hz is anticipated, 
f. Spectral Shape of Signal and Noise. Based on past experience, 
the signal and noise spectra are expected to peak in the range 
.5-to-3 Hz. 

Now that the assumptions about the data have been established, the 
choice of a suitable time-delay estimation approach can be made. The 
data statistical characteristics required for the Generalized Cross- 
Correlation method have been assumed in a, b, c, and d. Thus, the GCC 
method is applicable in this passive detection problem. As was shown in 
section 2.2, the basic cross-correlation has an undesirable spreading in 
its peak, which can be improved by a frequency weighting. The weighting 
proposed by Roth leads to a weighted correlation given oy 


V 2 (t) a s (T - T d> 


aG s s (f) 


j2irfT 


df 


Si % 


(f) 


+ G (f) 
n^ni! 


For the assumption made in f, about a similar spectral shape for sig- 
nal and noise [i.e., G (f) a const. G ss (f)], this correlation 
becomes 


X 1 X 2 


(t) 5 6 (t-t ) 


Thus, the Roth processor is appropriate for the assumed data, since its 
output is a delta function with a magnitude equal to one at the correct 
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delay. Implementing the Roth processor in the frequency domain, how- 
ever, has the disadvantages discussed in section 2.3. On the other 
hand, the advantages discussed in section 2.6 make the choice of imple- 
menting this processor in the time domain, using the IMS adaptive algo- 
rithm, evident. The choice of the TDE approach has been made. The 
question now is whether the adaptive filter discussed in section 2.5 is 
totally suitable for this application. This is the subject of the next 

r 

section. 

3.3 Filter Design 

The filter structure discussed by F. Reed, P. Feintuch, and N. 
Bershad [16], and presented in section 2.5, requires prior lead-lag 
information about the monitored signals in order to successfully esti- 
mate the delay. The delay estimate can be made only if the leading 
signal is fed as the primary input to the filter. This can be seen by 
examining Figure 2.3 where the input vector 

X T (n) = [x (n) x(n-l) x(n-2) ••• x(n-N+l)] 

component of the second sensor output, d(n) must lag x(n) so that 

the delay estimate can be made. 

In acoustic detection problems, the lead-lag information is 

unknown. This fact results in the necessity of applying the algorithm 

twice to the same data, swapping channels the second time, in order to 

get the lead-lag information as well as the delay estimate. This 

difficulty can be solved by introducing an — point delay in the 

2 
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secondary filter input as illustrated in Figure 3.1. Introducing this 
delay in d(n) leaves enough leading as well as lagging information in 
X(n) to be compared with d(n) to eventually produce the estimate. 

, In this implementation, the filter inputs will have a slightly 
different form than before. With discrete time signals given as in 
Chapter 2, 


Xi (n) = Si (n) + mi(n) 

X2(n) = % (n) + m2 (n) = a Si (n-D) + m2(n) 
the filter inputs are now given by 

x (n) = xi (n) = Si (n) + mi (n) (3.1a) 

and d(n) =x 2 (n-N') =a Si(n-N'-D) + m 2 (n-N') (3.1b) 

where 0< n< record length, |D| < N', N* = integer (— ) , and N is an 

2 

odd positive integer equal to the filter length. If, as assumed 
earlier, it^ (n) and n^n) are uncorrelated, having rr^n) delayed by 
N' points have no bearing on the analysis. All the analysis, presented 
in Chapter 2, of the weight vector is still applicable here except that 
D must be replaced by D+N 1 everywhere in the derivations. Equation 
(2.31) becomes in this case 
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W lms (K) = 0 sine {ir [K - (D+N)]} (3.2) 

As can be seen from (3.2), the weight vector now peaks at K = D+N'. So 
when 0=0, the weights peak at N' (the middle weight), and the peak 
is shifted right or left with positive or negative values of the delay, 
respectively. 

In section 2.4, it was shown that convergence of the weight vector 

f 

to the Weiner-Hopf solution can be obtained only if the stability con- 
stant K $ is bounded by equation (2.25), that is, 

-1 

Input power 

However, maintaining a constant K $ through the adaptation process can 
lead to instability in situations where a rapid change in the power 
level of the signal may occur. For this reason, K $ is updated in this 
design with every iteration of the algorithm. An estimate of tne input 
power of the signal, based on the signal components in the tapped delay 
line, is made with the arrival of each data sample. The value )K $ j is 
then chosen to be less than the inverse of the power estimate by a fact- 
or to be experimentally determined. Updating K $ in this fashion en- 
sured a fast, stable convergence of the weight vector regardless of the 
power level of the input. 

The final design issue is the convergence criteria. Due to the 
fact that our sampled data records may not contain any coherent signals 


-1 


< K < 0 


N 

Z 

i=l 


E[x?] 
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all the time, and due to the fact that the shape of the weight vector is 
unpredictable in these situations, a convergence criterion must be de- 
rived to distinguish between coherent and noncoherent input data before 
attempting to estimate the delay. This criterion is proposed as 
follows. Let us monitor the peak of the weight vector at every iter- 
ation as the data record is processed through the algorithm. By storing 
the abscissa value (corresponding to the discrete time delay) of the 
peak with every iteration, a data sequence is produced that can be plot- 
ted to produce a TDE vs. time plot. A variance criterion can then be 
applied to this TDE vs. time plot to determine whether a convergence 
has occurred or not. The data sequence representing tne abscissa values 
of the peaks is as long as the data record itself. Thus, for short-term 
stationary data, the variance criterion should be applied to a sliding 
window, of length less than the data records, through the sequence in- 
stead of computing the variance for the whole sequence at once. When- 
ever the variance criterion is met, the delay estimate is made based on 
the abscissa values contained in the window where the convergence has 
occurred. This provides a better estimate of the delay, compared to 
taking the abscissa value of the peak in the last iteration, especially 
in noisy environments. 

A detailed description of our implementation of the adaptive LMSTDE 
algorithm is presented next in Figure 3.2, Figure 3.3, and Figure 3.4, 
by means of flowcharts. In Figure 3.2, it is assumed tnat the sampled 
data records have been stored in arrays x and d of length L, where 
L is the data record length. W is an N point array, where N is 
the filter length, storing the weight vector components. The abscissa 
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Figure 3.3. Variance computation. 
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Figure 3.4. Delay estimate computation. 
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values of tne peaks will be stored in an L point array P. VR and 
VT denotes the variance and the variance threshold, respectively. The 
estimate of the stability constant, K , is computed with every iter- 
ation as 

K. = iL (3.3) 

Cu 

t 

where C is a predetermined factor and 

N-l N-l 

m = 2 E[x?(n)] = N x E [x 2 (n)] = I X 2 (n-i) 

i=0 i=0 

The peak detection block in Figure 3.2 is a simple search routine, wnile 
the variance computation is described in Figure 3.3 and tne delay compu- 
tation in Figure 3.4. The variance computation is made, as described 
earlier, over a sliding window of length K, where K < L. The equa- 
tion used for the variance estimate is 

i K-l , K-l 

VR = i e P 2 (i ) - [1 £ P(i)]2 (3.4) 

K 1*0 K 1*0 

Finally, the delay estimate is made by simply averaging the delay values 
contained in the window that yielded the minimum variance level. 

By inspecting the algorithm flowcharts, it is easily seen that the 
computational complexity is proportional to L x N. In practice L is 
the minimum number of iterations necessary for the convergence of the 
weight vector to the Wiener-Hopf solution. Therefore, for a specific 
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statistical type of data, L is a constant and the complexity becomes 
of order N. This design, accommodating positive and negative delays, 
implies a lower-bound on N of twice the maximum propagation delay 
between two sensors. Compared to the filter discussed in section 2.5, 
the filter length now is twice what it was before. With the complexity 
being of order N, it appears as if no speed improvement has been 
established over the previous design, where two passes were needed to 
obtain the delay estimate. Actual implementations, however, proved 
about IQ% time saving. The time saving was mainly obtained oy avoiding 
the overhead computation involved in every pass as can oe seen in the 
flowcharts. 

3.4 Specifications and Hardware 

This section presents and justifies the specific parameter values 
in the algorithm implementation as well as in the hardware system used. 
These specifications are summarized in Table 3.1. They can be divided 
into three categories. The first category includes parameter values 
derived from the assumption made in the previous section about the sig- 
nal frequency band. The second consists of parameters that were experi- 
mentally determined, while the ones related to hardware make up the 
third category. 

The starting point in the first set of parameters is tne frequency 
band tnat was specified as 1 to 20 Hz. An immediate result of this band 
is the sampling rate. A sampling rate of 40 Hz, satisfying the Nyquist 
rate, will not be adequate, the reason being that the signal has to be 
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Table 3.1. Parameter Specifications 


PARAMETER 

NOTATION 

VALUE 

Frequency Band 

— 

2 Hz - 20 Hz 

Sample Rate 

f s 

51.2 Hz 

Array Foot 

l 

800 ft. 

A/D 

— 

12 bits 

Quant. Level 

— 

2.4 mV/bit 

File Length 

M 

1024 points 

Record Length 

L 

128 points 

Filter Length 

N 

77 points 

Scale of Feedback 
Estimate 

C 

10 

Variance Threshold 

VT 

1 

Variance Window 

K 
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low-pass filtered at 20 Hz with a non-ideal roll-off filter. To elimi- 
nate any alias effects produced by attenuated frequencies aoove 20 Hz, a 
sampling rate of 51.2 Hz was chosen. The exact value of 51.2 Hz is 
desirable because it produces a rounded frequency resolution of .2 Hz on 
a 256-point FFT spectrum. Given the sampling rate, the array foot 
length (distance between two sensors) becomes a compromise between 
computational speed and bearing resolution. The shorter the array foot, 
the less the (maximum) delay, and the lower the resolution. The filter 
length, and hence the speed is directly proportional to the maximum 
propagation delay as was seen in the previous section. The bearing 
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resolution is related to the array foot as 

A9 - 90--COS-1 x speed 0f S0Und 1 (3.5) 

array foot length 

where A0 is the bearing resolution and At is the sampling period. 

To develop insight for the range of values involved, we can compute the 
filter length and bearing resolution for three different values of the 
array foot. (At = 1/51.2 sec, sound speed = 1080 ft/sec, filter 
length =■ 2 0 max +l) 

Array foot 2000ft. 800ft. 100ft. 

Filter length 200^ 77 # 10 o 

Bearing resolution .6° 1.5° 12° 

The first set of values implies a good resolution at the expense of slow 
speed while the third set implies a high speed computation at the ex- 
pense of poor bearing resolution. The compromise is obtained in the 
middle set. Note that the chosen distance between sensors is not likely 
to violate our assumption of uncorrelated noise across the array. This 
is true since infrasound noise produced mostly by local wind fluctu- 
ations and small eddies are not expected to propagate over an 800 foot 
distance. 

Experimental ly determined parameters are the record length (number 
of iterations), scale of stability constant, variance window, and vari- 
ance threshold. The first two parameters, and the second two, are re- 
lated. The record length, L, should be the minimum number of iter- 
ations that leads to convergence in the weight vector to the Wiener-Hopf 
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solution. The rate of convergence is controlled by the stability con- 
stant, K , as was seen in equation (2.21) in section 2.4. Keeping one 
of the two constant and varying the other, the system performance was 
observed and the values were chosen as 

L = 128 iterations 

C = scale of = 10 

The system performance, as L and C vary, is simulated in section 
4.2. The variance window, K, and the variance threshold, VT, were 
used in the convergence criteria to determine whether the data is coher- 
ent or not. Ideal response would have a zero variance across a 128- 
point (record length) window, but obviously that is too tight a criteri- 
on in practice. Based on the response of the algorithm to some experi- 
mental data gathered, including space-shuttle launch noise, moderate 
values of 48 points and 1 were chosen for K and VT, respectively. 

It is worth mentioning at this point, that although the algorithm 
computation is Dased on 128 iterations, the data files sampled are 1024 
points each. This file length is necessary for subsequent spectral 
analysis. 

A block diagram of the system hardware and most of its specifica- 
tions is presented in Figure 3.5. The type of sensors used in the array 
is the GLOBE 100 C Capacitor Microphone by GUS Manufacturing Inc. The 
high sensitivity (.2 volt/dyne/cm 2 ) and the low frequency response (.1 to 
(.1 to 500 Hz) on these Microphones make them very adequate for this 
application. Their dynamic range is 74 dB, .004 to 20 volts peak-to-peak 




Figure 3.5. System block diagram. 
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peak. This dynamic range implies an electrical noise amplitude of 
(upto) 4 m volts p-to-p. To avoid any electrical noise in the sampled 
data, the analog-to-digital converter should have a quantization level 
of at least 2 mV/bit when no amplification is used. But since a large 
dynamic range A/D is obviously desirable, a 12-bit, 2.4 mV/bit A/D con- 
verter was chosen. An analog input range of -5 V to +5V corresponds to a 
digital output range of 0 to 4095. Therefore, as long as no amplifi- 
cation is used on a microphone's output, the data are free from electri- 
cal noise. The specifications on the remaining hardware blocks are 
totally independent. The array configuration was intended to provide 
symmetry between the microphones so that delay-to-angle conversion can 
be made easier. Due to the experimental nature of tnis project and due 
to the lack of knowledge about the sound pressure level of infrasounds 
generated by meteorological events, the ampl ifier/attenuator block was 
necessary. Even though a signal frequency band of 1 to 20 Hz is antici- 
pated, the input is high-pass filtered at 2 Hz in order to eliminate the 
strong wind noise in the region of .1 to 2 Hz. Finally, the choice of 
the APPLE computer, as the system processing unit, was made to meet the 
economical, portable system requirement. 

3.5 Summary 

The design process started by evaluating the anticipated signal and 
noise. The main assumption made was that of similar spectral shape for 
signal and noise over a frequency band of l-to-20-Hz. This assumption 
led to the choice of the Roth processor which, in this case, produced a 
snarp peak at the correct delay in the weighted correlation function. 
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Examining the facts of section 2.6, proved that implementing the Roth 

processor in the time domain (rather than frequency domain as in the GCC 

method) using the LMS adaptive filter is advantageous. The specific 

N 

filter design for this application was then presented. An — -point 

delay in the secondary filter input was necessary to accommodate for 
positive and negative delays. Updating the stability constant, K , 
with every iteration was made for fast staole convergence of tne weight 
vector when the power level of the signal may change rapidly. To dis- 
criminate against non-coherent input signals in the automatic detection 
mode of the algorithm, a convergence criterion was devised. This crite- 
rion was based on the variance level of the delay estimate with time. 

The design phase ended by schematically presenting the exact implementa- 
tion of the design and specifying the software, as well as the hardware, 
parameters values. 



Chapter 4 


SIMULATION AND RESULTS . 

4.1 Introduction 

Aided by the theory of TDE techniques reviewed in Chapter 2, the 
system design has been described in the previous chapter. This chapter 
presents an evaluation of the foregoing design. The response of the TDE 
algorithm, as presented in section 3.3 is evaluated in section 2. Its 
response for deterministic data sequences is simulated as the different 
parameters are individually varied. Its effectiveness in tracking time 
delays between real, long-range, infrasounds is also illustrated. Sec- 
tion 3 describes the complete software package, with the TDE algorithm 
being its central processing element, used in the actual field operation 
of the system. Its routine operation and general results over the 
period of June to September 1984, are also discussed. The system is 
evaluated with respect to man-made and weather-related infrasounds. The 
man-made infrasounds are discussed in section 4, with an illustrative 
example using signals from jet-powered aircrafts, while the weather- 
related infrasounds are the subject of section 5. 

4.2 System Behavior 

This section is dedicated to the study of the system response for 
deterministic input data as well as real infrasound waves generated oy a 
far-field source. Deterministic data are used to study the system 
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Denavior as the different parameters, with effect on the response, are 
varied one at a time. The type of data used in the simulation and tne 
response for a nominal set of parameters are discussed in sub-section 1. 
The set of parameters can be divided into two groups. The first group 
consists of the adaptive filter parameters that were experimentally 
determined in the design. This group is the subject of subsection 2. 

The second group, containing the signal parameters (signal characteris- 
tics), is discussed in subsection 3. Finally the response for real 
signals consisting of the space shuttle launch noise is illustrated in 
the fourth subsection. 

4.2.1 Nominal Response 

The data used in this simulation are white data sequences generated 
on an Apple II computer using a random number generator function. White 
data are used for both signal and noise. Even though data with decaying 
spectra seem more appropriate for the simulation, white sequences were 
chosen for the ease in which they are generated. As explained in sec- 
tion 2.5, both types of data are expected to produce tne same response 
since the weight vector, for signal and noise having the same spectral 
shape, is a sine function regardless of what the spectral shape is. A 
256-point sequence and its corresponding auto-power spectrum, represent- 
ing the type of white sequences used, are plotted in Figure 4.1. The 
auto-power estimate is obtained by a 256-point FFT, and represented by 
the first 129 components since the last 127 are a mirror image of the 
first. The spectrum is expected to be completely white, in contrast to 
the one in the figure, as the number of overaged 256-point periodograms 
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approaches infinity. 

The nominal response of the algorithm for positive and negative 
delays is presented in Figure 4.2 and Figure 4.3, respectively. The 
positive delay in Figure 4.2 was generated in the input sequences as 
follows: 


x(n) 

= Ri (n) 

+ R 2 (n) 

(4.1a) 

d(n) 

= R]_ (n+4) 

+ Mn) 

(4.1b) 


where x(n) is the primary input to the filter, d(n) is the secondary 
input, R x ( n ) , R 2 (n) and R 3 {n) are uncorrelated random sequences. The 
embedded delay is four samples while the other parameters are set to 
their nominal values as specified in the figures. The negative delay 
was obtained by simply swapping the adaptive filter inputs, i.e.: 

x(n) = Ri(n+4) + R 3 (n) 
d(n) = Ri (n) + R 2 (n) 

Tne fact that x(n) now lags d(n) by four samples leads to a negative 
delay in the output of the filter as illustrated in Figure 4.3. 

The first plot in both figures was obtained by monitoring the 
abscissa value of the peak in the weight vector at every iteration of 
the algorithm. The second plot is the final set of weights. The 
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4.2. Nominal filter response for D = +4 samples. (C=l 
N = 17, L = 128, SNR = 1). (a) TDE vs. number of 

iterations, (b) Final weight vector. 
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N=17, L=128, SNR=1) . (a) TDE vs. number of iterations 
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capability of tracking positive and negative delays is evident from the 
figures. The "steady-state" delay has shifted from +4 to -4 in the TOE 
vs. time plot, while the peak in the weight vector has shifted right and 
then left by 4 units as expected. These plots are to be compared with 
the response obtained in the next two subsections where the various 
parameters are varied. 

/ 

4.2.2 Filter Parameters 

The filter parameters include the filter length (N), stability 
constant scale factor (C), and the record length (L). In the design 
process, a lower bound on N was set while C and L were experiment- 
ally determined. This simulation verifies the necessity of tne lower 
bound on N and explores the basis on which C and L were deter- 
mined. The input sequences used in this subsection are the ones previ- 
ously given in section (4.1). 

The lower bound on tne filter length (number of weights) was de- 
rived in section 3.3 as (2xD +1). Twice the maximum delay was neces- 

' max 

sary since the filter accommodates positive and negative delays. A 
number of weights less than the lower bound will obviously give no esti- 
mate. A very large number of weights, on the other hand, will not 
improve the response since the sine function (the shape of the weight 
vector) approaches zero as the number of weights is increased. This 
behavior is exemplified for D=4 in Figure 4.4 with N=7 and in Figure 
4.5 with N=33. The remaining parameters are set to their nominal 
values in both figures. The response in both cases is to be compared 
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Figure 4.4. N=7. Delay can not be estimated with N<2xD 

(D=4 , C=10, L=128, SNR=1) . (a) TDE vs. ma 

number of iterations, (b) Final weight vector 
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with the nominal response in Figure 4.2 where N=17. 

The effect of varying C on the stability of the system is demon- 
strated in the next two figures. It was seen in section 3.3, the bounds 
on are given by: 


< K < 0 

N-l S 

2 X 2 ( T ) 

' 1=0 

The value used in this implementation is given as 

K - = — 

5 N-l Cp 

C 2 X 2 ( 1) 
i=0 

where C is a positive scale factor that has to be greater than, one to 
satisfy the stability condition. The output of the filter is plotted in 
Figure 4.6 and Figure 4.7 for C=.5 and C=200, respectively. Having 
C less than one led to an unstable system as predicted by Widrow [11], 
while increasing its value behind ten hardly affected the response. 

The record length (number of iterations) was chosen as 128. The 
purpose of this simulation is to show that a 128-iteration process is 
sufficient for the convergence of the weight vector to the Wiener-Hopf 
solution. This is illustrated in the next two figures in comparison 
with Figure 4.2. Figure 4.8 illustrates the ambiguity in the filter 
output caused by too small a number of iterations. The expected 
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Figure 4.8. L=32. Ambiguity in filter output due to short record 

length. (D=4, N=17, SNR=1). (a) TDE vs. number of 

iterations, (b) Final weight vector. 
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response, after a satisfactory convergence, was obtained for L=128 in 
Figure 4.2. A larger number of iterations produced slightly improved 
results, though not appreciable, for the additional computation. This 
aspect is shown in Figure 4.9 for L=512 points. 

4.2.3 Signal Parameters 

This subsection illustrates the system response as the cnaracteris- 
tics of the filter inputs are changed. The characteristics considered 
here are signal-to-noise power ratio (SNR), signal frequency bandwidth, 
and time varying delays. 

The signal to noise ratio used in the nominal response of Figure 
4.2 was SNR=1. The response is evaluated in the next two figures for 
SNR =.5 and SNR=2. The input sequences were generated as 

x(n) = SxRi (n) + R 2 (n) 

d(n) = SxRi(n+4) + R 3 (n) 

where x(n) is the primary input to the filter and d(n) is the sec- 
ondary one. Rj ( n ) , R^n), and R 3 (n) are uncorrelated random sequences 
with the same amplitude. The scale factor S equals .5 in Figure 4.10 
and 2 in Figure 4.11. In all three cases, the response of the system 
was satisfactory after 128 iterations. However, as can be seen from 
the figures, the lower the SNR the slower the response (convergence). 

The next simulation demonstrates the dependency of the main lobe 
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Figure 4.9. 1=512. Improved response over L=128 but not necessary. 

(0=4, C=10, N=17, SNR=1) . (a) TDE vs. number of iterations 

(b) Final weight vector. 
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Figure 4.10. SNR=l/2, Slow Convergence. (0=4, C=10, N=17, L=i28). 

(a) TDE vs. number of iterations, (b) Final weight vector. 
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Figure 4.11. SNR=2. Faster convergence. (D=4, C=10, N=17, L=128) . 

(a) TDE vs. number of iterations, (b) Final weight vector 
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width of the weight vector on the signal frequency band. The expression 
for the weight vector for an ideal band-limited signal and white addi- 
tive noise was derived by Ahmed and Carter [17] as 

W lm$ (K) = 46w b sinc[o> b (K-D)] cos [^(K-d)] 

where 3 is a real constant, to, = , to = 4 L Tr^ iL . and w o 

d 2 c 2 

are the signal frequency limits, i.e. G ss (to) = 0 for <o < to 0 and u 

> u)^ , where Ul > co 0 . It was shown in section 2.5 that the width of 
the main lobe of the weight vector is given by — , which becomes 

O) (0 1 

C 1 

for ojq = 0. The input sequences were generated as 


X ( n) = S x Rj (n) + Mn) 


d(n) = S x Ri(n+4) + R 3 (n) 


where x(")> d(n), Rz(n), and R 3 (n) are the same as before, S is a 
scale factor used to obtain a SNR=1, and R (n) is an ideally filtered 

sequence specified by its amplitude spectrim as 


Ri (to ) 


_ r 1 for 0 < <0 <0)1 

‘ 0 for oai < o> < it 


In Figure 4.2, the frequency band of Ri(n) was given by too = 0 and 
<o 1 = tt , i.e. white sequence. The main lobe in that case was 2 
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resolution ( — = 2) wide. For the simulation in Figure 4.12, the 

Oil 

bandwidth of Ri(n) is limited to — , i.e. 

4 

IT IT 

G D D (o) ) = 0 for — < a) < tt -*■ u) 0 = 0 and . 

K i K i 4 4 

As seen in the figure the main lobe width now is = 8 resolutions. 

it 

t ______ 

4 

Finally, u> 0 = 0 and = — were used in Figure 4.13. Producing a 

8 

main lobe that is 16 resolutions wide. 

It was mentioned in section 2.6 tnat one of the advantages of the 
LMSTDE algorithm over other methods is that it has the aoility to tracK 
moving sources. The routine is tested here for this ability by linearly 
varying the delay between the input signals in a noisy environment. The 
input sequences were generated as 

x(n) = Ri(n) + R 2 (n) 
d(n) = Ri(n+D(n)) + R 3 ( n) 

where Ri (n), R 2 (n) and R 3 (n) are uncorrelated random sequences as 

before, and D(n) is a time varying delay. The signal-to-noise ratio 
used here is 2. In Figure 4.14, with D(0) = 0, D(n) was increased by 
one every 64 points in 896-point input sequences. The solid line 
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.13. Wi= 1 . Main lobe width of 16 resolutions. (D=14, 010 
N=33, L=128, SNR=1). (a) TDE vs. number of iterations, 

(b) Final weight vector. 
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Figure 4.14 Filter output for time varying delay (1 point per 64 points). 
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represents the actual delay embedded, while the dashed line represents 
the instantaneous estimate of the delay as the abscissa value of the 
peak in the weight vector during 896 iterations. Figure 4.15 differs 
from Figure 4.14 in that D(n) is varied at the rate of 1 point every 
128 points instead of 64. The adaptive filter tracked more closely the 
actual delay in the second case. This is to be expected, since the 
weights were given more 'time to converge. In both cases however, D(n) 
represents a relatively fast varying delay simulating a fast moving 
source. Considering this fact, one can conclude that the system 
response is acceptable and that the time lag, predicted by Feintuch, 
Bershad, and Reed [18], between the estimate and the actual delay is not 
as severe as it appears. 

4.2.4 Shuttle Launch Noise 

Long-range infrasound signals produced by a large rocket launch 
have been recorded and studied by several investigators [22]. These 
signals are similar in their statistical characteristics to weather- 
related infrasounds. Thus, they can be used to simulate the anticipated 
infrasounds generated by a far-field meteorological event. In this 
subsection, the adaptive filter algorithm is applied to the space 
shuttle launch signals and the overall response is evaluated. 

Space Shuttle VIII was launched from Cape Canaveral, Florida, at 
2:31 a.m. EST on 30 August 1983. About 52 minutes later the launch 
signal arrived at Wallops Island, Virginia, and was sampled on an array 
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Figure 4.15. Filter output for time varying delay. (1 point per 128 points). 
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of four infrasonic microphones at the rate of 14 Hz. The array configu- 
ration is illustrated in Figure 4.16. Time histories and auto-power 
spectrums, of one channel before, during, and after the event, are pre- 
sented in Figure 4.17 and Figure 4.18, respectively. These two figures 
represent a verification of the event and show the background of the 
signal that can be considered as the uncorrelated noise on the different 
microphones. The adaptive filter algorithm was applied to three pairs 
of channels during a 9.14 sec. (128 points) time frame. The output of 
the microphones during this time frame is plotted in Figure 4.19. The 
adaptive filter output for C2 x C3 (channel 2 cross channel 3), C3 x C4, 
and C4 x C2 is shown in Figure 4.20, Figure 4.21, and Figure 4.22, 
respectively. Making the estimate of the delays based on a 48-point, 
zero variance, window in the TDE vs. time plots yields the following 
results: 


C2 leads C3 by .455 sec. (+7 resolutions) 

C3 lags C4 by .780 sec. (-12 resolutions) 

C4 leads C2 by .325 sec. (+5 resolutions) 

Based on these delays and on the array configuration, a source direction 
of about 210° from north can easily be obtained using trigonometric 
identities. This estimate of the Dearing is to be compared with tne 
true direction of 207° as illustrated in Figure 4.23. The error in the 
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Figure 4.17. 


Time histories of STS-8 launch signal (in equivalent 
arbitrary units) on one channel, (a) before, (b) during 
and (c) after the arrival of launch signal. 
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Figure 4.19. Time histories of STS-8 launch signal (in equivalent 

arbitrary units) on three microphones, (a) channel 2, (b) 
channel 3, (c) channel 4. 
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Figure 4.20. Filter response for launch signal, C2 x C3, ( 
vs. time, (b) Final weight vector. 
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Figure 4.21. Filter response for launch signal, C3 x C4, (a) TDE 
vs. time, (b) Final weight vector. 
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estimate is within the bearing resolution. In this case, the bearing 

r i 

resolution (for f $ = 14 Hz and array foot of 650 ft) is about 6.8°. 

4.3 System Operation and Results 
The previous section demonstrated the effectiveness of the LMS 
adaptive filter algorithm in estimating the delay between the outputs of 

two sensors. The algorithm, as discussed in section 3.3, can now be 

/ 

used as a subroutine in a larger software package for a complete, actual 
field operation of the system described in section 3.4. This section 
presents a description of this software package, the mode in which it is 
operated, and the way the results are evaluated. 

The software package is depicted in Figure 4.24. It consists of an 
automatic detection routine, as well as a bearing computation, of infra- 
sonic sound waves. The whole process is automatic in the sense that no 
user intervention is needed after the program is started and the param- 
eters are set. However, this automation is limited by the data storage 
capacity of the computer. The process is started by simultaneously 
sampling three 1024-point files from the three sensors* in the array. 

The data files are stored in temporary buffers until they are analyzed 
by the LMSTDE algorithm. The algorithm, as described in section 3.1, is 
then applied to each pair of data files at a time, considering only the 
first 128 points as the filter data record. If the convergence test, 
using the variance criterion, fails for any pair of sensor outputs, the 
process will be stopped and a new set of data files is sampled. When 


*0nly three out of the four sensors in the array are at this stage. 
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all three pairs of files show good coherent data through acceptable 
variance levels, the data files are transferred to a floppy disk for 
permanent storage. 

As the detection phase ends, the spectral content of the signal, 
the bearing, and the nature of the source are then evaluated. The 
analysis begins by testing the delays, obtained from the LMSTDE algo- 
rithm, for consistency assuming plane wave propagation and 0° elevation. 
If the consistency test is not passed the delays are printed for later 
analysis; otherwise three angles, corresponding to the three delays, are 
evaluated using simple trigonometric relations. Tne delays and angles 
are then printed for the record, and the whole process is started over 
again. 

The consistency test discriminates between plane waves and cylin- 
drical, or spherical propagating waves. The later category may occur if 
the source is in the near field of the array (within one mile from the 
center of the array). Since the intention of this research is the early 
detection of far-field events, this category is neglected by the soft- 
ware. The delays obtained in such cases are only printed for direct 
analysis by the user. The trigonometric relations used to evaluate the 
bearing in case of a plane wave were formulated for horizontally incom- 
ing waves only. In this case, the three computed angles will be the 
same while they will differ from each other for a plane wave propagating 
with a non-zero elevation. The difference among the angles can be 
easily related to the value of the incident angle. The software does 
not produce this calculation. For further discussion see Appendix A. 
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The software has been applied since early June 1984 at NASA Langley 
Research Center. An average of ten data files have been stored daily, 
about half of which indicate plane wave propagation with printed delays 
and angles. This printed data is then compared with the weather re- 
ports, obtained from the meteorology center at NASA, on a daily basis. 
Finally, the data files stored on floppy disks will undergo a comprehen- 
sive spectral analysis, njainly for signature identification and power 
law verification. The software used for the spectral analysis is, or 
has stemmed from, DAISE (Digital Analysis of Infra-Sound Experiments) 
[23], that was developed previous to the course of this research. 

Daily evaluation of the data proved that most files stored by the 
system can be grouped into four categories, with their typical sources, 
as follows: 



Man-made 

Weather-rel ated 

Near-field* 

Heavy machinery 

Microbursts 

Far-field 

Military aircrafts 

Severe storms 


*Within one mile from the center of the array. 

Man-made infrasounds were to be expected since the microphone array is 
located in Langley Research Center where a wide range of activities take 
place on a regular basis. These include operating all sorts of machin- 
ery, wind tunnels, etc. Future plans for near-field (both categories) 
signal analysis exist but they are not of interest at this point. 
Infrasounds produced by far-field events are discussed in the next two 
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sections. The man-made category is discussed in section 4.4 while the 
weather-related category is discussed in section 4.5. 

4.4 Man-Made Signals 

This section discusses the man-made signals with long-range propa- 
gation detected by the automatic detection routine described in the 
previous section. This category of infrasounds is of interest to this 
study for two reasons. First, it verifies the capability of the system 
to detect infrasonic signals and to locate their sources. Second, it 
illustrates the necessity of more sophisticated signal processing tech- 
niques to distinguish between these events and tne meteorological ones. 

Among others, the non-weather related signals detected include 
single-component infrasounds (for example 7.2 Hz, 10 Hz, 16.8 Hz), sonic 
booms generated by supersonic vehicles, and infrasounds from jet powered 
aircraft. Single component infrasounds have been detected daily. The 
system estimate of their directions places their sources in the main 
research center of Langley. Thus, it is believed that they are mostly 
related to the different activities in the Center. Few sonic booms were 
detected during the summer of 1984. They were mainly associated with 
flights of supersonic military aircraft such as the Air Force F- 15 which 
are stationed nearby at Langley Air Force Base. The high altitude of 
these sonic booms was indicated by elevation angles of 25° to 40° evalu- 
ated from the data. The third type of man-made signal was found in 
several files over a period of about three months. The common charac- 
teristics between tnese infrasounds include a sound pressure level of 
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about 90 dB, frequency band of about 1 to 6.5 Hz, and a source direction 
ranging from 130° to 150°, relative to true north. The consistency of 
these values between files, the relatively strong amplitude, and the 
directions found all have led to the belief that these infrasounds are 
generated during the take-off of military jet powered aircraft from 
Langley Air Force Base. The location of the runway at the base relative 
to the microphone array -is shown in Figure 4.25. 

One of the jet data files is considered here for illustration. The 
file was stored on August 28, 1984, by the automatic detection routine. 
The parameters of the software and the hardware were set to their nomi- 
nal values as specified in Table 3.1. The signal on one channel before 
and during the event is shown in Figure 4.26 in the time domain and in 
Figure 4.27 in terms of its auto-power spectrum. The signal on micro- 
phone 1, 2, and 3 during a different time frame, of the same file, is 
shown in Figure 4.28. To demonstrate the convergence of the algorithm 
for this type of data, the filter response was recomputed for Cl x C2 in 
Figure 4.29, C2 x C3 in Figure 4.30, and C3 x Cl in Figure 4.31. This 
response was recomputed by feeding alternative data points of the time 
histories to the adaptive filter. This effectively reduced the sampling 
rate to 25.6 Hz (51.2 Hz originally). A 25.6 Hz sample rate is adequate 
in this case since the signal has a narrow frequency band as can be seen 
in Figure 4.27. Keeping this rate in mind, examining the filter 
response yields the following results 


Cl lags C2 by 0.742 sec. (-19 resolutions) 




Figure 4.25. Location of the runway in Langley Air Force Base relative 
to the microphone array. 
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Figure 4.2b. Time histories (in equivalent arbitrary units) of infra- 
sounds generated by jet-powered aircraft, (a) before, and 
(b) during the event. 
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4.27. Auto-power spectrums (in arbitrary units) of infrasounds 
generated by jet-powered aircraft, (a) before, and (b) 
during the event. 
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Figure 4.28. Time histories (in equivalent arbitrary units) of infra- 
sounds generated by jet-powered aircraft, sampled on (a) 
channel 1, (b) channel 2, and (c) channel 3. 
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Figure 4.29. Filter response for jet data, Cl x C2, (aj TDE vs 
time, (b) Final weight vector. 
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C2 leads C3 by 0.625 sec. (16 resolutions) 

C3 leads Cl by 0.117 sec. (3 resolutions) 

Based on these delays, a source direction of about 145° can be found. 
The estimate of the bearing places the source roughly on the middle of 
the runway as seen in Figure 4.25. This result among others obtained 
from similar files presents a reasonable assurance that the system is 
indeed capable of detecting and locating the sources of infrasounds. 

4.5 Weather-Related Data 

As was mentioned in Chapter one, the far-field meteorological 
events of interest in this research are clear-air turbulence and severe 
storms. There has not been clear evidence of clear-air turbulence 
detection on the system since it has been in operation. It is also 
believed that only a few storms were recorded. The uncertainty govern- 
ing their detection is high in all cases except for one severe storm. 
This storm was the hurricane "DIANA." It was relatively close (200 
miles) and stationary for at least one day. The data files stored by 
the system during that day, and their evaluation, present enougn neces- 
sary and reasonably sufficient conditions to identify the data with the 
hurricane. An overall evaluation of the system, including a discussion 
of the various factors that might have contributed to the lack of mete- 
orological data, is discussed in the conclusion chapter. This section 
will concentrate on the hurricane and the issues surrounding its detec- 


tion. 
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A tropical storm developed into a hurricane over the coastal water 
of North Carolina during the period of September 10 to September 12, 
1984. It reached its full strength by September 12 with recorded wind 
velocity up to 110 mph. The location of its "eye" on that day was rela- 
tively stationary about 200 miles away from NASA/Langley Research Center 
at latitude 33.8 and longitude 77.3. During the afternoon and evening 
of the same day, the system recorded eleven files, four of which pointed 
toward the site of the hurricane with a calculated sound pressure level 
of 85 dB. The first three files, were stored within approximately a 40 
minute time frame. This raises the possibility of having all three 
files associated with the same event. Examination of the spectral shape 
and contents of the signal on one channel through all four files sup- 
ports this possibility. For reference, information in the signal sam- 
pled on one microphone is presented in Figure 4.32 as time history, in 
Figure 4.33 as auto-power spectrum, and in Figure 4.34 as log-amplitude 
vs. log-frequency auto-power plots.* Compared to other data files re- 
corded, these data present three, somehow unique, features. First, the 
signal has a wide frequency band. This can not be seen on the linear 
scale of Figure 4.33, because of the difference in the dynamic range of 
the components on both ends of the spectrum. But it is clear on the 
log-scale of Figure 4.34 where the band plotted is 2-to-16 Hz. Second, 
the roll-off of the spectrum is consistent with the prediction of 

-3 5 

Meecham and Ford [21]. They predicted a roll-off proportional to co , 

*Software used to obtain these plots was developed at NASA/Langley 
Research Center, by Dr. Allan J. Zuckerwar and David Katzoff. 



Discrete time (in seconds) 

Time histories (in equivalent arbitrary units) of hurricane signal sampled on one 
channel at different times, (a) File one, (b) File two, (c) F,le three, (d) File four 


Figure 4.33. 


Discrete frequency (in Hz) 

Auto-power spectrums (in arbitrary units) of hurricane signal sampled on one channel at 
different times, (a) File one, (b) File two, (c) File three, (d) File four. 


log-ampl itude 




Figure 4.34. 


Auto-power spectrums on a log scale (log amplitude vs. log frequency, in equivalent 
arbitrary units) of hurricane signal sampled one channel at different times, (a) 
File one, (b) File two, (c) File three, (d) File four. 
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while the ones computed in Figure 4.34 range between w * and w * . 

Third, the frequency band as well as the roll-off were consistent 
through all the files even though the fourth file was sampled hours 
after the first three. 

To illustrate the consistency between the microphones outputs, the 
signal on all three microphones during the same time period was plotted 
in Figure 4.35. Its auta-power spectrum is plotted in Figure 4.36 on a 
linear scale and in Figure 4.37 on a log-scale.* The adaptive filter 
algorithm was reapplied to the data in Figure 4.35 to demonstrate the 
type of convergence obtained for this signal. The response is shown in 
Figure 4.38, Figure 4.39, and Figure 4.40 for Cl x C2, C2 x C3, and C3 x 
Cl, respectively. As can be seen in the figures, the TDE vs. time plots 
reflect the noisy (actually windy) environment in which the signal trav- 
eled and was then sampled. The directions evaluated from the different 
four data files were, with respect to the bearing resolution, within an 
acceptable tolerance from each other. A mean direction of 187°, from 
north, was obtained. This bearing is to be compared with the actual 
value of 193° as illustrated in Figure 4.41. 

The foregoing analysis and results clearly constitute the necessary 
conditions for true detection of the hurricane. The sufficiency, how- 
ever, is lacking a higher detection rate, on behalf of the system, over 
a longer period of time. This issue, among others, will be discussed in 
Chapter 5. 

*Software used to obtain these plots was developed at NASA/Langley 
Research Center, by Dr. Allan J. Zuckerwar and David Katzoff. 
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Figure 4.35. Time histories (in equivalent arbitrary units) of 

hurricane signal simultaneously sampled at (a) channel 
one, (b) channel two, and (c) channel three. 
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Figure 4.39. Filter response for hurricane data, C2 x C3, (a) TDE vs. time 
(b) Final weight vector. 
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Figure 4.40. Filter response for hurricane data, C3 x Cl, (a) TDE vs 
time, (b) Final weight vector. 



Figure 4.41. Hurricane site relative to the microphone array in 
NASA Langley Research Center. 









4.6 Summary 


The response of the basic LMSTUE algorithm was evaluated for deter- 
ministic data sequences. The simulation justified the specific choice 
of the filter parameters. The effect of varying the signal character- 
istics on the response was demonstrated. A direct proportional ity be- 
tween the SNR and the rate of convergence of the weight vector was ob- 
served. The dependency of the main lobe width of the weight vector on 
the frequency band of the signal, as predicted by Ahmed and Carter [17], 
was verified. The ability of the algorithm to track time-varying delays 
was also illustrated in the simulation. A steady state time-lag between 
the actual varying delay and the estimate of about 64 samples was ob- 
served for SNR=2. The algorithm was then tested using Space Shuttle 
VIII launch noise. The response was satisfactory in the sense that the 
launch site was pin-pointed within the bearing resolution. Hie specific 
software used in the field on a daily basis was described. It consisted 
of an automatic detection and bearing evaluation routine. Tne operation 
of the system over the period of June to September 1984 was evaluated 
with respect to the rate and type of signals detected. Data evaluation 
proved that man-made signals dominated the data files stored by the 
system. A sample of detected signals from jet-powered aircraft was 
presented. This sample verified the capability of the system to detect 
infrasounds and locate their sources. The system detection of weather- 
related infrasounds was not conclusive. Only one severe storm was de- 
tected with reasonable assurance. The assurance of the detection was 
based on the consistency of the signal characteristics and bearings 
through several files sampled on the day of the storm. 



CHAPTER 5 


CONCLUSION 
5.1 Remarks 

The goal of this research was the design of a system for detection 
and bearing evaluation of infrasounds. Based on the results presented 
in Chapter 4, it is believed that this goal has been achieved. The 
adaptive LMSTDE algorithm employed in the design was evaluated in sec- 
tion 4.1. Simulation results were in total agreement with theoretical 
expectations. The response of the algorithm to the shuttle noise clear- 
ly demonstrated its ability to estimate delays in cases of long-range 
infrasonic propogation. The automated- system performance, in the actual 
field, proved effective through detection and bearing evaluation of man- 
made signals. This class of signals consisted of several types of 
infrasounds including that of jet powered aircraft, whose signature is 
similar to what is anticipated in meteorological events. Acceptable 
performance of the system, as an infrasonic detector and bearing evalu- 
ator, was also achieved for large-scale weather events such as the 
hurricane "DIANA." 

This work, however, is faced with the difficulty of evaluating the 
system as a turbulence detector for two reasons, both of which are 
beyond the designer's control. First, the data obtained in this case is 
poorly defined. That is, it is not repeatable, nor is it unambiguously 
verifiable. The second reason is the lack of a large ntmber of weather- 


116 



117 


related data files needed for the evaluation. As pointed out earlier, 
few meteorological events were detected on the system. Only the hurri- 
cane was identified with reasonable assurance. Some of the factors that 
might have contributed to the low detection rate of events, or of files 
within the same event, are listed below. 

1. Low sound pressure level of weather-related infrasounds. A 
minimum level of 54 dB, derived from the system calibration, is 
required at the site of the array for the signal to be detect- 
ed. The fact that the hurricane, despite being a powerful 
event, produced only 85 dB, raises the question of how big and 
how close the event should be in order to be detected. This 
question is left open for acoustical investigators. 

2. Noisy environment. The presence of strong man-made signals at 
the array location makes it difficult for the system to identi- 
fy small-scale weather events. 

3. Wind effects. Small wind fluctuations at the array can be 
considered uncorrelated noise and usually do not affect tne 
system response. However, this is not the case in windy con- 
ditions. Strong wind may correlate at two or more microphones, 
destroying the coherence property of the turbulence signal. 

4. Multiple delay. In situations of scattered thunderstorms or 
large frontal systems, the signal sampled at the array is the 
superposition of infrasounds from more than one event. The 
multiple delay pattern embedded in this type of signal creates 
an ambiguity problem that can not be solved on the existing 
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system. 

r> 

5. Short-term stationary signals. Although the lifetime of a 
meteorological event is relatively long (several hours in some 
cases), the period over which the turbulence, and hence the 
infrasonic emission, exists is virtually unknown. Therefore, 
the event may be easily missed in the case of turbulence that 
is stationary ov,er a short time relative to the computational 
speed of the system; i.e. 20 seconds worth of data are sampled 
every 7 minutes (average computational time), so infrasonic 
emission for less than 7 minutes can be missed. 

5.2 Future ulork 

More work is needed to improve the system performance for turbu- 
lence identification purposes. One work area is the investigation of 
signal processing techniques that are appropriate for elimination of 
non-weather data and removal of wind effects from the signal before 
applying the adaptive LMSTDE algorithm. Oice enough knowledge is 
obtained about the infrasonic signature of turubulence, pattern recog- 
nition methods can be used to discriminate against non-weather data. 
Removal of wind effects may require the development of a wind-measure- 
ment system that is synchronized with the infrasonic detection scheme. 
The purpose of this system would be to provide the wind information 
needed for compensation in the system input. 

Another future research consideration is. the study and modification 
of the LMSTDE algorithm for supressing the side lobes in the weight 
vector and processing of multiple delay signals. Both issues are 



119 


important for avoiding any ambiguity in the weight vector peak that may 
be caused by the characteristics of the input signal. 

Future work also needs to consider implementing the LMSTDE algo- 
rithm in hardware using dedicated signal processing chips. Hardware 
implementation would dramatically improve the computational speed, and 
hence, the detection rate, of the system in the case of short-term 
stationary signals. 
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APPENDIX A 


This appendix represents the derivation of a relationship between 
the source bearing and source elevation, and the delays obtained for an 
arbitrary sound wave. A general relation is first derived for one pair 
of microphones. The specific relations for the array used in this work 

t 

are then discussed. 

Figure A.l illustrates the elevation effects on the delay between 
two microphones, where 

9 is the source direction relative to the line joining the 
microphones 

<t> is the source elevation angle 

D' is the delay due to the elevated-source signal 

D is the projection of D 1 on the horizontal plane (delay for <j>=0) 

V is the sound velocity 

z is the distance between microphones. 

Inspection of the figure yields the following two relations. 


and 


cos $ 


D'xV = IV 
DxV D 


(A.l) 


cos 9 


DxV 
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(A. 2) 
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substituting for D in (A.l) gives 


\i n 1 

cos <j> = _ x (A. 3) 

4 COS 0 

Equation (A. 3) is a transcendental equation that relates the direction 
and elevation of the source to the delay between two microphones. Note 
that D' is given, as thre output of the LMSTDE algorithm, while <j> and 
e are to be found. 

The array configuration, and the values of 0 ( 0 l5 9 2 > 0 3 ) for the 

three pair of microphones, are shown in Figure A. 2 (horizontal plane). 

It! 

Let , C> 2 , and D 3 be the delays obtained from microphones (3, 4), 
(3,2), and (3,1), respectively, due to a signal propogating with <t>° 
elevation. Applying equation (A. 3) to the aifferent pairs yields 


cos 


V 
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iL 


COS 01 


COS <f> 


I 



% COS 02 


I 



a cos 03 


(A. 4a) 


(A. 4b) 


(A. 4c) 


01 , © 2 , and 03 represent the source direction relative to the 
microphones axis as shown in the figure. If 6 represents the source 

direction relative to true North, then 0i , 02 , and 03 can be expressed 
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Figure A. 2. Relative angles. 
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in terms of g by linear relationships, i.e., in figure A. 2, the re- 
lations are 

9i * * - 180° 

0 2 =300°. 3 (A. 5) 

0 3 = 240°- B 

/ 

1 I I 

Hence, having , D 2 > and D 3 estimated by the system, the problem of 
finding <|> and 0 becomes that of solving two equations (out of the 
three in eq. (A. 4)) with two unknowns. 

Note that the relations of (A. 5) hold for a certain range of 3 
only. Thus, when solving for the elevation and the exact bearing of the 
source in a given situation, these relations have to be evaluated with 
the bearing values provided by the system. 



